import sys, os
from os.path import expanduser as home
import numpy as np
import scipy
from scipy.io import wavfile
from matplotlib import pyplot as plt

sys.path = ['..'] + sys.path # this is where the features module is stored
from features import logfbank, fbank

%pylab inline

WAV_DIR = home("~/Dropbox/Shared/Orthoptera Sound App/species_recordings")

Populating the interactive namespace from numpy and matplotlib

Utilities

class NotAWaveFileError(Exception):
    def __init__(self, filename):
        self.filename = filename
    def __str__(self):
        return "%s is not a WAVE file" % self.filename

def spectr_wave(filename, signal, fs, melfft, NFFT):
    '''
    Small waveform, large spectrogram
    '''
    
    fig = plt.figure(figsize=(15,5), dpi=300)
    fig.subplots_adjust(hspace=0.2,)
    
    ax1 = plt.subplot2grid((6,2), (4,0), rowspan=2) # waveform
    ax2 = plt.subplot2grid((6,2), (0,0), rowspan=4) # spectrogram
    ax3 = plt.subplot2grid((6,2), (0,1), rowspan=6) # FFT of spectrum
    
    ### PLOT 1 - Waveform      
    t = np.linspace(0, len(signal)/float(fs), len(signal))
    
    # if necessary, downscale the waveform
    ax = ax1
    d_signal, d_time = downscale(signal, t)

    ax.tick_params(labelleft='off')
    ax.autoscale(tight=True)
    ax.set_ylabel("Amplitude", labelpad=22)
    ax.plot(d_time, d_signal, rasterized=True)
    ax.set_xlabel("Time (s)")

    ### PLOT 2 - Spectrogram
    ax = ax2
    specgram = ax.specgram(signal, Fs = fs, scale_by_freq=True, rasterized=True)
    
    ax.autoscale(tight=True)           # no space between plot an axes
    #ax.get_xaxis().set_visible(False) # remove the x tick for top plot
    yticks = np.arange(0,fs/2,5000)    # |
    ax.set_yticks(yticks)              # |
    ax.set_yticklabels(yticks/1000)    # |
    ax.set_ylabel("Freq (kHz)")        # |> change Hz to kHz
    ax.tick_params(labelbottom='off') # labels along the bottom edge are off
    
    ax = ax3
    ax.pcolormesh(melfft)
    ax.set_xlim( (0,NFFT/2) )
    
    
    fig.suptitle(filename)

    # save
    plt.show()
    plt.close()
        
    # print "finished %s" % name
    
def stats(data, fs):
    fs = float(fs)
    nsecs = len(data)/fs
    lengths.append(nsecs)
    fss.append(fs)


def open_audio(filename):
    '''
    open audio file with scipy or audiolab, depending on what is available.
    Return data, fs, enc
    '''
    
    if not filename.lower().endswith(".wav"):
        raise NotAWaveFileError(filename)
    try:
        fs, data = wavfile.read(filename)
    except ValueError:
        fs, sampwidth, data = wavio.readwav(filename)
        print "\t%s is a %dbit audio file" % (os.path.basename(filename), 8*sampwidth)
        
    enc = None
    #data, fs, enc = audiolab.wavread(sound) # same with audiolab
        
    # if the audio sample is stereo, take only one channel. May not work as
    # desired if the two channels are considerably different. 
    if len(data.shape) > 1:
        data = data[:,1]
    
    return (data, fs, enc)

def downscale(data, time):
    data2 = []
    t2 = []
    q = 100
    r = 2*q
    for i in range(0,len(data),r):
        data2.append(max(data[i:i+r]))
        data2.append(min(data[i:i+r]))
        try:
            t2.append(time[i+int(0.25*r)])
            t2.append(time[i+int(0.75*r)])
        except:
            #print i+int(0.75*r)
            pass
    #plt.gca().get_xaxis().set_visible(False)
    data2 = data2[:len(t2)]
    return data2, t2

NFFT = 128

for soundfile in os.listdir(WAV_DIR):
    if not soundfile.lower().endswith(".wav"): continue
    signal, fs, enc = open_audio(os.path.join(WAV_DIR, soundfile))
    
    mel_feat = fbank(signal, nfilt=40, winstep=0.0232, winlen=0.0232)[0]
    melfft = abs(scipy.fft(mel_feat.T, n=NFFT))
    
    spectr_wave(soundfile, signal, fs, melfft, NFFT)